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Abstract. We generalize the technique of Lada et al. (1994) to map dust column density through a molecular cloud 



(Nice) to an optimized multi-band technique (Nicer) that can be applied to any multi-band survey of molecular 
clouds. We present a first application to a ~ 625 deg 2 subset of the Two Micron All Sky Survey (2MASS) data 
and show that when compared to Nice, the optimized Nicer technique (i) achieves the same extinction peak 
values, (ii) improves the noise variance of the map by a factor of 2 and (iii) is able to reach 3 a dust extinction 
measurements as low as Av — 0.5 magnitudes, better than or equivalent to classical optical star count techniques 
and below the threshold column density for the formation of CO, the brightest H2 tracer in radio-spectroscopy 
techniques. The application of the Nicer techniques to near-infrared data obtained with a 8 meter-class telescope 
with a state-of-the-art NIR camera, such as the VLT-ISAAC combination, will be able to achieve dynamic ranges 
from below 10 21 protons cm -2 to over 10 23 protons cm -2 (Ay in the range [0.3, 60]) and spatial resolutions < 10", 
making use of a single and straightforward column density tracer, extinction by interstellar dust. 
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1. Introduction 



in the range [0.5, 3] magnitudes, i.e., from column densities 
10 21 and 



between ~ 1 



10 21 protons cm 2 (Rossanc 



the Un verse. Allhough Ihey are ubiquitous ill galaxies like 198g . ' gee u & Goodman 19991 for a multi-technique 



the Milky Way, we sLill do not have a basic understanding — approach) 
of how they relate to the more diffuse interstellar medium 
(ISM), of how some of them form stars and planets, or how 
they vanish. We do know today that they are the coldest 
objects known (T ~ 10 K), which unfortunately implies 
that they are also one of the hardest to observe and study. 



Dark cta-orfe remain urn of Lhe lcasL uudersLood objects in 197g |. |Di ckman l 97 8|; |Mattila 1986|; fcregorio Hctem^I 



Almost everything we know today about molecular 
cloud structure has been derived from radio observations 
of molecular tracers (e.g . , CO, CS, NH 3 ) of the unde- 
tecta ble H 2 QLada 199^ ; |Blitz fc Williams 1999] [Myers 



Be fore the discovery of emission from molecules in dark 1999| ) and in more recent years through dust thermal 



clouds ( Wilson et al. 1970] , after which these clouds be- 
came known as molecular clouds), the general technique 
used to study the density structure of these objects relied 
on statistical analysis of numbers of dust extincted stars in 
the background of a cloud. This approach, introduced by 
Wolf (1923), is known as the star count method and was 
later refined and intensively applied by Bok (1937, 1956| ). 
In this technique, the main mass component of a molecular 
cloud, H2, is traced by the dust in the cloud, an assump- 
tion later proved to be valid by observations. The main 
feat of the star count method is its straightforward essence 
and its sensitivity to low column density regions with Ay 



emission that can be detected by radio continuum tech- 
niques (Andre et al. 2000: Johnstone et al. 200C). Radio 
spectroscopy and continuum techniques are able to probe 
deeper into molecular clouds, well past column densities 
of 10 22 protons cm -2 , an order of magnitude improve- 
ment over classical star count techniques. Moreover, radio- 
spectroscopy offers an unique view of the dynamical struc- 
ture of molecular clouds (e.g., Rosolowsky et al. 1999), 



something that none of the absorption or emission contin- 
uum techniques can, evidently, assess. However, the inter- 
pretation of results derived fr om radio observat i ons is not 
always straight forward (e.g. Alves et al. 1999 ; |Chandlci 



fc Richer 2000] ). Several poorly constrained effects inher- 
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ent to these techniques (e.g., deviations from local ther- 
modynamic equilibrium, opacity variations, chemical cvo- 
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lution, small-scale structure, depletion of molecules, un- 
known emissivity properties of the dust, unknown dust 
temperature) make the construction of an unambiguous 
picture of the physical structure of these objects a very dif- 
ficult task. The price of being able to probe deeper into a 
molecular cloud with radio line and continuum techniques 
is a complicated data analysis, not free from ambiguities. 

It has long been recognized that infrared color excess 
can be used to obtain a reliable estimate of the extinc- 



tion th rough a molecular cloud (see, e.g., |Bok 1977| ; [Jones 
et al l|98(f [Frcrking et al. 1982[ [Jones ct al. 1984| ; |CasaIi 
198fj[)~TIic deployment of sensitive, large format infrared 
array cameras on large telescopes, however, has allowed 
the direct measurement of line-of-sight dust extinction to- 
ward thousands of individual background stars observed 
through a molecular cloud. Such measurements are free 
from the complications that plague molecular-line or dust 
emission data and enable detailed maps of cloud density 
structure to be constructed. In all fairness, all density 
tracer methods rely on the assumption of an universal 
gas-to-dust ratio, a fair assumption on theoretical grounds 
that still remains to be demonstrated for cold dense molec- 
ular clo.nHs Larla et al (1994) pionpprerl a t.errmiqnr-; the 



Near-IMrarerl Color Eycpss (Nice) technique; for mea- 



suring and mapping the distribution of dust through a 
molecular cloud using data obtained in large-scale, multi- 
wavelength, infrared imaging surveys. This method com- 
bines measurements of near-infrared color excess to di- 
rectly measure extinctions (as opposed to the statistical 
measurement in the star count method) and map the dust 
column density distribution through a cloud. Moreover, 
the measurements can be made at significantly higher an- 
gular resolutions and substantially greater optical depths 
than previously thought possible. The efficacy of this tech- 
nique was demonstrated with the study of the dark cloud 
L 977 ( [Alves et al. 1998|), IC 5 146 ( |Lada et al. 1999| ), and 
Barnard 68 ( Alves et al. 2001 ) where through straightfor- 
ward analysis of nearly 7000 infrared sources background 
to these clouds produced detailed maps of the extinction 
to optical depths and spatial resolution an order of mag- 
nitude higher than previously possible (Ay in the range 
[1,40] magnitudes, spatial resolution down to 10 arcsec). 



1.1. The large view on molecular clouds 

Nearby galactic molecular cloud complexes represent our 
best chance to understand molecular cloud structure. 
These complexes have relatively large sizes (from a few 
to tens of parsec) and are hence fairly extended, typically 
stretching several to tens of degrees across the sky. A large 
view of entire molecular cloud complexes, with sufficient 
resolution and dynamic range, is the only way to put these 
clouds in their Galactic context and address the questions 
on their origin and fate. Census of entire molecular cloud 



dynamical ranges of 10 21 protons cm -2 to below or about 
10 22 protons cm -2 , but suffer from the caveats inherent to 
the respective used techniques (see above) . There is a clear 
need for a new large scale mapping technique of molecular 
cloud complexes, that makes use of a more reliable column 
density tracer, and is able to offer an extended dynamical 
range. 

In this paper we generalize Lada's Nice technique 
to an optimized multi-band technique, the Near-Infrared 
Color Excess Revisited technique (Nicer). Although in- 
s pired by the Two M icron All Sky Survey (2MASS) 
(Klcinmann ct al. 1994), and by the possibility of all sky 
dust extinction mapping, this technique can be applied to 
any multi-band survey of molecular clouds [e.g., DENIS 
( |Epchtein et al. 1997| ), SDSS ( |York et al. 2000| ), or any 
combination of multi- wavelength catalogs]. We present a 
first application to 2MASS data, and show that Nicer 
improves the noise variance of a map by a factor of 2 
when compared to Nice, and is able to reach 3 o dust 
extinction measurements of Ay ~ 0.5 magnitudes, bet- 
ter than or equivalent to optical star count techniques 
and below the threshold column density for the for ma- 
tion of CO under typical t emperatures and de nsities ( van 
Dishoeck fc Black 198^ ; |van Dishocck 1992] ). This lat- 
ter point is highly relevant as CO is the brightest H2 
tracer in radio-spectroscopy techniques. The application 
of Nicer to near-infrared (NIR) data obtained with an 
8 meter-class telescope outfitted with a state-of-the-art 
NIR camera, e.g. the VLT-ISAAC combination, will be 
able to achieve column density dynamic ranges from below 
10 21 protons cm -2 to over 10 23 protons cm~ 2 (Ay in the 
range [0.3,60] magnitudes), on spatial resolutions < 10", 
making use of a single and straightforward column density 
tracer, extinction by dust. 

The paper is organized as follows: In Sect. 2 we review 
the Nice technique and discuss its limitations, in Sect. 3 
we give a detailed presentation of Nicer, in Sect. 4 we 
apply both techniques to the same 2MASS data subset, 
and present the conclusions in Sect. 5. 



2. Near-infrared color excess 

The difference between the observed and the intrinsic color 
of a background star provides information on the opti- 
cal depth of the cloud. For example, fo r the near-infrared 
and a ssuming a normal reddening law ( Rieke fc Lebofsky 
1985), we have 



(H - K) obs = {H - K) tT + 0.063 • A x 



(1) 



complexes are rare (CO line emission: [Dame ct al. 2001 
Fukui ct al. 1999| ; IRAS opacity maps: [Wood fc Myers 



1995| ; optical star counts: |Cambresy 1999] ), have typical 



where H and K represent the H and K-band (centered 
at 1.65/mi and 2.2/mi, respectively) the superscripts "obs" 
and "tr" denote observed and true (intrinsic) quantities, 
respectively. The visual extinction Ay is often used as a 
measure of dust column density and will be the main tar- 
get of our investigations. This qua ntity is relevant because, 
as suggested by several studie s ( Lilley 1955 ; Jenkins fe 
[Savage 1974 ; Bohlin ct al. 1978 ), it is directly related to the 
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projected mass density of the cloud. In particular, it has 
been shown that most clouds have the same gas-to-dust 
ratio, so that N(H2) — Ayl.9-10 21 protons cm~ 2 mag -1 
or, equivalently for the projected two-dimensional density, 

K = Ay ■ 15 Mq pc~ 2 . 

Equation (|]) would not be particularly useful if the 
intrinsic star colors (which are generally unknown) would 
span a large range of values. In reality, infrared colors of 
stars have a relatively small scatter (e.g., a{H — K) tr ~ 
0.08), so that it is sensible to take the average of Eq. (|l|), 
obtaining 



-4, 



15.87- 



{(H - K) ohs ) - {(H ~ K) tT ) 



(2) 



It is not difficult to estimate ((H — K) tT ) using, for ex- 
ample, a control field where the extinction is negligible. 
Hence, in principle, we can use Eq. (§) with ((H-K) ohs ) 
replaced by the observed color of each background star and 
obtain an unbiased estimate for the cloud column density 
in the direction of the star. In practice, this method suf- 
fers from two main limitations: (i) The use of the observed 
color of each star in lieu of ((iJ — K) ohs ) introduces a 
source of noise in the estimate of Ay; (ii) It may not be 
easy to establish if a given star is a background object, 
i.e. is farther than the cloud. Traditionally, both limita- 
tions are dealt with by using two simple techniques. The 
noise on the measured column density, which is ultimately 
introduced by the intrinsic scatter on the colors of stars, 
is reduced by averaging angularly close objects (which, 
hopefully, are subject to the same amount of extinction; 
see however below). Thus, normally the extinction is esti- 
mated using the equation 



Ax 



15.87 



N 



n=l 



(3) 



where the sum is performed on a set of N angularly close 
stars. This is the appr oach taken by the Nice method, 
originally described by Lada et al. (1994 ). Regarding the 
second problem, foreground stars can be recognized from 
their color, which is significantly bluer than angularly 
nearby stars (often the number of foreground stars is much 
smaller than the number of background stars) . 

It is not difficult to recognize that the two problems 
discussed above are intimately related. Suppose, for ex- 
ample, that we are investigating a cloud with low column 
density. Then, because of the intrinsic scatter of star colors 
and of photometric errors, we could take a foreground star 
as a background object just because its is observed redder 
than expected. The situation can be even more difficult 
for embedded stars. Clearly, a perfect discrimination be- 
tween foreground and background objects is feasible only 
in the unrealistic case where the column density can be 
measured without any error. This is not possible, but at 
least we can try to take care of both problems in a uniform 
way. This is the point of view adopted by our optimized 
method for extinction estimate. 



Fig. 1. JHK color-color of one hundred stars (sketch). 
The left plot shows the colors of stars which are not sub- 
ject to extinction. The ellipse encloses a 3cr confidence 
region; note that the ellipse is almost exactly vertical ori- 
ented, i.e. J — H correlates only weakly with H — K. The 
right plot shows the same stars as observed through a 
cloud with Ay = 5 magnitudes (a normal reddening law 
is assumed). 



3. The optimized method 

In this section we present Nicer, the Near-Infrared 
Color Excess method Revisited, an optimized, multi-band 
method for extinction estimation based on the near in- 
frared excess of background starlight. The method has 
been designed by considering separately the two steps 
needed to make an extinction map: 

— Local extinction estimate for each star; 

— Spatial smoothing of individual stars estimates. 



We describe these two steps in Sects. 3.1 and 3.2 



3.1. Local extinction 

As discussed in Sect. |[ the cloud column density Ay is 
normally measured by comparing the H — K color for stars 
observed through the cloud with the same color for stars 
for which no reddening is expected [see Eq. (||)]. Clearly, 
other possibilities are also viable. For instance, we could 
choose the J — H color and use the expression 



Ay = 9.35 



1 N 

n=l 



(4) 



Should we use Eq. (||) or (|J) to infer the cloud column 
density? In both cases, we expect some error on the esti- 
mate of Ay. The final error on Ay depends critically on 
three factors: 

1. The scatter of the intrinsic colors: Normally, stars 
present a larger scatter in J — H than in H — K (see 
Fig. |). 
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Fig. 2. For each star, Nicer obtains an estimate of the ex- 
tinction Ay using a maximum-likelihood technique. The 
star colors and the relative covariance (represented in the 
figure in the left panel by the small ellipse to the top right) 
are obtained from the observations. Note that correlation 
of colors is expected [see Eq. dip))]. The net effect of pho- 
tometric errors is a widening (solid ellipse) of the intrin- 
sic color scatter (dashed ellipse) [Eq. ([n]) and following- 
discussion]. To estimate the extinction to the star Nicer 
than finds the ellipse center along the reddening line (see 
figure to the right) that is closer to the observed colors of 
the star [Eq. @]. 



2. The photometric error on the individual bands: 
Typically errors on J are smaller than errors on K (ex- 
cept along lines-of-sight through dense clouds, where 
extinction makes stars fainter in J). 

3. The numerical coefficient used to convert color excess 
into Ay extinction: This coefficient is larger for Eq. (||) 
(15.87 vs. 9.35). 

Hence, on the one hand we are encouraged to use Eq. (||) 
because of the smaller scatter in colors; on the other hand 
we want to avoid this estimator because of the larger nu- 
merical coefficient and larger photometric errors. In con- 
clusion, this Heuristic discussion does not lead to a defini- 
tive answer regarding the choice between the two estima- 
tors (||) and (Q). This suggests that a more quantitative 
approach is needed to study the problem and decide which 
estimator should be used. Actually, the analysis we are 
going to carry out will provide a much more interesting 
result and will show that a clever use of all available in- 
frared bands produces significantly more accurate results 
than the individual estimators (||) or ([|) . [In the following 
we will consider the typical case of three bands J, H, and 
K available; The generalization to more bands is trivial.] 
With the three bands at our disposal we can make two 
independent colors, c\ = J — H and c 2 = H — K . For each 
Ci (with i = 1,2), we can write the relationship between 
the observed color c° bs and the true one c' r (i.e., the color 
that would be observed if no extinction were present) as 

c° bs - cf + hAy + a . (5) 



Here ki — Ei/Ay is the ratio between the color excess on 
the band i and the extinction on the V band; for the colors 
considered here we have ki = 1/9.35 and k 2 = 1/15.87. 
The extra term a above represents the noise on the colors, 
i.e. the result of photometric error on the estimate of J—H 
and H — K. Let us restrict the discussion to estimators of 
Ay which are linear in the observed colors c° bs , i.e. of the 
form 

A v = a + &icf s + b 2 cf s . (6) 

The coefficients a, b\, and b 2 need to be determined so as 
to satisfy two conditions (see Fig. ||): 

1. The estimator is unbiased, i.e. its expected value is the 
true extinction Ay. 

2. The estimator has minimum variance. 

The first condition is equivalent to the equations 

hh + b 2 k 2 = 1 , a + &i(ci r ) + 6 2 (c 2 r ) = . (7) 

The variance of the estimator Ay can be evaluated from 
the expression 

Ay - (Ay) = J2 b * " (<?)) + E b ^ ■ ( 8 ) 

i i 

As a result, we can immediately write 

Var(iy) = b * h 3 Covy (c tr ) + ^ b * b 3 Covy (e) . (9) 

The covariance matrices introduced in this equation rep- 
resent two different sources of noise. The first matrix, 
Covy(c tr ) is the intrinsic scatter of star colors, which 
clearly makes the determination of the extinction Ay un- 
certain (for example, if a star is redder than expected 
we will overestimate Ay). This matrix, obviously inde- 
pendent of the star considered, can easily be determined 
by using a control field, where the extinction can be ne- 
glected, and measuring the scatters of star colors. 

The second matrix, Covy(e), is related to photometric 
errors, and thus changes for each star. This matrix can be 
easily calculated provided that an estimate of the photo- 
metric errors for each star is available. In the typical case 
of the three bands J, H, and K considered here, we can 
write 

cov« w = (i+r» a . (10 ) 

Note that, while the photometric errors on different bands 
(namely, aj, an, and ok) can be taken to be uncorrelated, 
errors on colors are not. 

We can now minimize Eq. ([)]) with the requirement 
that Eqs. (Q) holds using Lagrange's multipliers. This way 
we reduce our problem to the solution (61, b 2 ) of the linear 
system 

Cn C 12 -fcA fh\ /0\ 

C21 C 22 -k 2 U 2 = , (11) 
-hi -hi / \\J 
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where the matrix CVj = CoVy(c tr ) + Covy(e) is the sum 
of the two covariances. The parameter A in the previous 
equation is Lagrange's multiplier and, for our purposes, 
its value is uninteresting. The solution of this system is 
then, in matrix notation, 



where is the weight for the n-th star, given by 



C" 1 k 
k -C- 1 k ' 

where b = (61,62) and k = (fci,/c 2 ) 
for the optimal estimator is thus 



(12) 



The final expression 



A v = b 1 [c° hs - (<$)] + b 2 [cf 



(13) 



with b given by Eq. (|12|) 
or equivalently Var(yly) 



Its variance is given by Eq. 
= b-Cb. 



A n ice feature of this technique is that it can be applied 



withoult significant modifications in the case where one of 

the bands is not observed. Suppose, for example, that the 
J band is not available. Then we can assume an arbitrary 
value for this band and set the corresponding error o~j to 
an extremely large value. If we "blindly" apply Eqs. ( |l2|) 
and (|Hi|), we obtain a matrix C _1 with only a single non- 
vanishing element at the position (2,2). In other words, 
the use of a large error on J automatically suppresses the 
use of this band in the evaluation of Ay . Interestingly, the 
same technique can be used if the missing band is H , which 
contributes to both colors c\ = J — H and c 2 = H — K: 
In this case the estimator will be composed only of the 
combination c\ + c 2 = J — K , thus avoiding the use of the 
missing H band. 

3.2. Spatial smoothing 

So far we have considered only a single star. In order to 
obtain a smooth extinction map with high signal-to-noise 
ratio, Nicer applies a spatial smoothing to angularly close 
stars. The smoothing is performed using three different 
techniques, described below in individual subsections. 

The choice of the smoothing technique is important 
for two main reasons: (i) The final map has a signal-to- 
noise ratio that strongly depends on the smoothing used; 
(ii) The smoothing is responsible for the selection of back- 
ground stars. As already pointed out, Nicer tries to deal 
with both problems at the same time. 



3.2.1. Weighted mean 

The simplest way to obtain a smooth mass map from in- 
dividual estimates for Ay is to use a weighted mean. More 
precisely, we use the values of Ay obtained for stars close 
to a direction 9 on the sky to estimate the extinction on 6. 
Calling 0( n ) the position on the sky of the n-th star, Ay 



w{e-ew 

Var(A 



(«)> 



(15) 



Note that, in contrast to Eq. (g), here the index n runs 
over all N observed stars in the field. The weight W(") is 
the combination of two different factors: (i) A spatial term, 
described by a weight function W; (ii) An error weight, 
proportional to l/Var(^V). The weight function is usu- 
ally chosen to be a Gaussian with appropriate width. We 
stress that the choice of the characteristic length of the 
weight function is a fundamental step in the reconstruc- 
tion process, as shown by a detailed statistical study of the 
statistical properties of the smoothed map. This has al- 



ready been do ne in a different context ( Lombardi fc Bcrtin 
1998| ; see also Lombardi fc Schneider 20011) , and thus we 



just state the main results obtained. The size of the weight 
function directly determines the correlation length of the 
final extinction map, or, equivalently, sets the scale of the 
smallest details we can hope to observe on the map (every 
feature smaller than the characteristic size of the weight 
function will be washed out). Moreover, the larger the 
width of the weighting function, the larger the effective 
number of stars used for each point on the map, and the 
smaller the error on the extinction. In other words, we 
can decide whether we want a high signal-to-noise ratio 
map with low resolution, or a more noisy map with higher 
resolution. A local estimate for the error a s is given by 



£Li(^ W ) a Var(4 B) ) 



(16) 



where Var(yly^) is the expected variance for the n-th star 



[see note after Eq. (13)]. Finally, we note that the choice 
( |l5"| ) for the coefficients of the stars maximizes the signal- 
to-noise ratio. 

So far we have assumed that all stars are background 
objects with respect to the cloud. In reality, some of them 
could be in the foreground and could thus contaminate 
our extinction estimate. If the number of foreground stars 
is not too large, we can ignore them and proceed as if all 
stars were in the background. Because of the linearity of 
the estimator (|l4|), we can correct for the dilution of the 
signal a posteriori by multiplying the map Ay (6) by the 
factor A7iVb ac k, where A^ack is the estimated number of 
background stars. Actually, this simple technique is often 
not very efficient and c an als o introduce biases for dense 
clouds (see below Sect. 3.2.4 ). For this reasons, it might 
be preferable to use other smoothing techniques described 
below. 



its estimated extinction, and Vax(A v 
estimated variance, we write 



Ay {9) 



the corresponding 3.2.2. Sigma-clipping 

Foreground stars do not contribute to the signal, still con- 
tribute to the noise. Hence, it is desirable to exclude in the 
(14) average ( |l4| ) those stars that are suspected to be in the 
foreground because of their low extinction. This can be 
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easily obtained by using a sigma-clipped or by a median 
estimator (described below in Sect. 3.2. 3| ) . 

In order to implement sigma-clipping, we make a first 
estimate of Ay using Eq. (|l4|) and calculate the expected 
variance from Eq. (^). Then we repeat the entire pro- 
cess using only those stars for which Ay^ does not dif- 
fer by more than a factor 3(7^ from the first estimate 

of A v evaluated at (n \ i.e. Ay (9^). This procedure 
quickly converges to a pair of values, the measured ex- 
tinction Ay (9) and the related error, which we take as 
fiducial values. 

In order to verify the goodness of our procedure, we 
can evaluate the observed variance on the data: 



Elt^] 2 



(17) 



Note that this variance is defined in a different way with 
respect to Eq. ( |l6|) . Assuming that the photometric errors 
on the various bands have been correctly evaluated, we ex- 
pect (7 ; 



> o ^ r , the two variances being equal only if all 



perform the spatial smoothing. In some conditions, sev- 
eral subtle problems, in fact, may make each technique 
described above inappropriate, noisy, or even biased. 

As already noted by Thoraval et al. (1997 ), the sim- 
ple or weighted mean (as described in Sect. 3.2.1 ) can be 
severely biased in high-extinction regions. In fact, because 
of extinction, the density of background stars decreases as 
the cloud column density increases, while the density of 
foreground stars is left unchanged. As a result, in the dense 
regions foreground stars represent a larger fraction of the 
local number of stars with respect to less dense regions. 
This way, we could under-estimate the column density in 
the central regions of clouds. Such a bias can generally be 
ignored for very nearby clouds. In other cases, the bias 
might be reduced using the median or the sigma-clipping 
methods. 

On the other hand, substructure on the scale of the 
weighting function can also be a problem. In case of 
the simple weighted mean, the smoothing procedure is 
straightforward to understand (at least at the simple level 
kept in the discussion here; see however ILombardi & 



stars a:-e subject to the same extinction, i.e. if there are no Schneider 2001): The final map obtained is expected to be 



foreground stars and the cloud is perfectly homogeneous. 
Actually, the difference between dv| and a ^ can be used 
to obtain precious information on the substructure of the 
molecular cloud (see, e.g. Juvela 199^). 



3.2.3. Weighted median 

Alternatively, we can use a weighted median in order to 
make the extinction map. This quantity is defined as the 
solution A(9) of the non-linear equation 



N 



J2 W {n) sgn(i(6>) - iW) = 



(18) 



where sgn is the sign function (sgna; = {+1,0,-1} de- 
pending on the sign of x). Since in general the condition 
( |l8| ) cannot be satisfied exactly, a linear interpolation is 
used (see Lombardi et al. 2000] for a similar example of 
application of a weighted median). 

The use of the median has several advantages with re- 
spect to the simple mean. Probably, the most interesting 
is the fact that the median is a robust estimator, since it 
is basically unaffected by outliers. This property is par- 
ticularly useful in our case, since it provides an efficient 
method to remove foreground objects from our map. On 
the other hand, the median has the significant disadvan- 
tage of having noise properties difficult to study. We also 
note that the implementation of the median is not as effi- 
cient as the mean, but this should be considered a minor 
problem. 

3.2.4. Discussion 

In the previous subsections we have described three simple 
smoothing techniques used by Nicer. An important point 
to note is that there is not really a single best method to 



the original, true extinction of the cloud convolved with 
the weight function used. In the case of sigma-clipping 
or median smoothing, things are much more complicated. 
Suppose, for example, that the cloud does not present any 
significant substructure on the scale of the weighting func- 
tion, except for small "holes" with very low absorption. In 
this case, we would identify stars observed through these 
holes as foreground objects, although they are in reality 
in the background. As a result, we would completely miss 
the holes, and thus overestimate the cloud column den- 
sity. Similarly, for a cloud with substructures consisting of 
dense globules, we would underestimate the column den- 
sity (note that both sigma-clipping and median are "sym- 
metric," in the sense that they discard outliers with both 
large or small Ay). 

In summary, each smoothing technique has advantages 
and disadvantages (this is the ultimate reason to include 
all of them in Nicer) . A key role is played by the fraction 
of foreground stars (and thus the distance of the cloud) , 
and by the substructure of the cloud. It is in any case 
always advisable to use all smoothing techniques and to 
check the consistency of the maps obtained. 

4. Nicer at work: The Orion region from 2MASS 
data 

The Two Micron All Sky Survey (2MASS) offers a unique 
opportunity to test our algorithm. Conceived about a 



decade ago (Kleinmann et al. 1994), the 2MASS project 
is an all-sky-survey obtained from two 1.3-m telescopes in 
the J (1.25 /im), H (1.65 //m), and K s band (2.17 /jm). 
Although the surveys has been completed, only about 47% 
of the sky is at present publicly available in the Second 
Incremental Release.^] 

1 See http://www.ipac.caltech.edu/2mass/. 
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We selected from the 2MASS archive a large area 
around the Orion and MonR2 nebulae, with galactic co- 
ordinates 

200° < I < 230° , -30° < b < -5° . (19) 

The Orion region is a good candidate for our study for 
several reasons. The density of background stars, although 
not particularly high, is enough to perform a good analy- 
sis. Moreover, the 2MASS second incremental release cov- 
ers this region very well, with only a small number of 
"holes" (except for the Nord-West part). Finally, there 
are already a number of studies in the literature which 
can be used to test our results (e.g., Carpenter 2000 uses 
the same data set to investigate the spatial distribution of 
young stars in these clouds and presents averaged J — K s 
color maps of this region). 

For each point source in this coordinate range we re- 
trieved the magnitude on the three bands (fields x_mag), 
the relative total errors (x_msigcom), and some additional 
flags used to classify the goodness and the reliability of 
the detection (cc_f lg and rd_f lg). Possible artifacts (ob- 
jects with cc_f lg 7^ 0) were excluded from the analysis. 
We then applied our dedicated pipeline, as described in 
the following sections. 



4.1. Preliminary analysis 

Given the area on the sky and the number of stars, the 
pipeline "suggested" to perform the analysis using a grid 
of 820 x 700 points, with scale 2.4' per pixel, and with 
Gaussian smoothing characterized by FWHM = 4.8'. By 
default, the pipeline uses a FWHM of 2 pixels, and about 
15 "effective" stars per point (i.e., about 15 stars con- 
tribute to the signal for each point). We then made a pre- 
liminary analysis using standard values for the average 
colors and color scatters, thus obtaining a first extinction 
map for the region. The map obtained was used only as a 
first check for the basic parameters chosen, to identify a 
control field, and to obtain the photometric properties of 
stars in our sample. In particular, inside the boundaries 
, we found a region in the South-East of the map which 
does not show any significant extinction and which thus 
has been selected as control field (see Fig. |]). We then 
used the color of stars in this control field to set the pho- 
tometric parameters needed in our algorithm, namely the 
average colors of stars (i.e., (J — H) and (H — K)) and 
the scatter in colors (or, more precisely, the covariance 
matrix, which describes also correlations between colors). 
This step was performed in our pipeline by flagging stars 
inside the control region (which can be defined using the 
standard S AO image format for region files) and by making 
a statistical analysis on the photometric data of flagged 
stars. Note that here only the 2MASS photometric data 
were used, and not the preliminary map (this map was 
used only to select the control field). As a result, an inac- 
curate choice for the initial constants of Eq. (Q) to make 
the preliminary map had no influence on the statistical 
parameters obtained. 



Fig. 3. Top: Color-color diagram for stars in the control 
field, made from the colors of 6,000 stars. The ellipse rep- 
resents a 3cr confidence region. Bottom: Color-color dia- 
gram for stars in the central region of the cloud. The red- 
dening due to the dust in the cloud significantly stretches 
the locus of stars in the plot. 



We observe that the use of a large survey such as 
2MASS has the significant advantage with respect to nor- 
mal observations of letting us study in detail the halos of 
molecular clouds. As a result, the choice for the control 
field is very robust, in the sense that even regions with 
extremely low extinction can be avoided. We also stress 
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Fig. 4. The total region mapped using Nicer. The cuts used in Ay emphasize the faint, diffuse halos around the main 
cores. The low noise of this map allows us to detect a Ay = 0.5 extinction with 3 a confidence level. The higher noise 
observed in the southern part is due to the smaller star density. 



that this is possible because of the availability of data on 
large areas. Fig. || shows the color dispersion of stars in 
the control field, together with the 3a elliptical confidence 
region. 



4.2. Final map 

The final extinction map, shown in Fig. was performed 
using the method discussed in this paper. Stars in the 
control field was used to set the intrinsic colors and color 
dispersions of (i.e., to evalua te t he quantities (cj) and 



Cov ij (c tr ) described in Sect. For each star in our 

region, we evaluated the column density using Eqs. ( fl2|) 
and (|l^) . Finally, column densities for close stars was av- 
eraged out using three different schemes, namely simple 
(weighted) mean, sigma clipping, and median. 



Figures f| and || show the result obtained using a 
Gaussian smoothing characterized by FVHM = 5' with 
er-clipping (other smoothing schemes gives similar results, 
and differences cannot generally be seen by eye). Note that 
use of our optimized method, combined with the quality of 
the 2MASS data, allows us to detect at an unprecedented 
level of detail the faint, extended halos on which the main 
cores are embedded. At the smoothing size used, we reach 
the excellent sensitivity Ay = 0.13 at la level. 

Since the focus of this paper is mainly on the opti- 
mized reconstruction method, we will defer a discussion 
of the results obtained on this and on similar cloud com- 
plexes to a future paper. However, we want to point out 
a couple of peculiar features of the map obtained. First, 
we detect a single region with significantly bluer color (see 
Fig. |§|, right). This region corresponds to the open clus- 
ter NGC 2204, whose blue stars contaminate our map thus 
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Fig. 5. Zoom of the central region of the Nicer extinction map. Cuts in this map emphasize the faint, diffuse halos 
around the main cores. Contours are displayed at Ay = {2, 4, 8} magnitudes. Note that the maximum value obtained 
for Ay is 17.7 magnitudes. 



simulating a "negative" extinction. Note also the cometary 
structures observed south to the Orion cloud (Fig. [|, left). 
In these condensations we observe extinction Ay > 6 on 
scales close to our smoothing length. Clearly, it would be 
worthwhile to study these regions at higher resolution. 

In order to better appreciate the advantages obtained 
using the optimized technique described in this paper, we 
also carried out the analysis on the same field using only 
the weighted average on the H — K color for stars with 
relatively low photometric errors (we required errors on 
H and K to be smaller than 0.15 magnitudes). Figure [7] 
shows the result obtained with this standard method, and 
should be compared with Fig. S. 

— We first note that, apart from the different noise level, 
the two maps basically give the same values for the 
column density. This proves that the new method de- 
scribed in this paper is reliable and unbiased. 



— The simple H — K map is significantly noisier than the 
one obtained using the optimized method. Note that 
the resolution, set by the scale of the Gaussian smooth- 
ing, is the same in both maps. A quantitative analysis 
shows that with the use of the optimized method we 
gain a factor of two on the variance of the map. In 
other words, the new technique is able to reach the 
same signal-to-noise ratio of the standard method us- 
ing only half of the stars. 

— Because of the reduced noise level, in Fig. [|we are able 
to distinguish regions with extremely low column den- 
sity {Ay ~ 0.5 magnitudes), which would be otherwise 
undetectable. 

We note that the exact gain obtained with the optimized 
method strongly depends on the individual intrinsic scat- 
ter of stars on each color. This last point depends on sev- 
eral factors, such as the depth of observations, the galac- 
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Fig. 6. The noise estimate [cf. Eq. for the central region shown in Fig. ||. The average noise ranges from 

about 0.15 to the top of the figure to about 0.25 to the bottom (where the density of background stars is smaller 
because of the increase in |6|); contours are plotted at cr^ = {0.25,0.35}. In the white regions the error diverges 
because of the lack of data. Increase in noise is observed in high absorption regions (because of the decrease in density 
of stars), and around bright stars in the top of the figure (which have been masked by the 2MASS pipeline). The 
white dot on the right is due to bright star zeta Ori (V = 1.79 mag). Except for these features, the overall noise level 
is rather uniform. 



tic latitude of the field, and the colors used. On the other 
hand, given the excellent results obtained with the 2MASS 
archive, we are confident that Nicer can produce extinc- 
tion maps with significantly higher signal-to-noise ratio 
for any infrared data. Finally, we stress that although 
the method has been presented here for the J, H, and 
K bands, it can actually be applied to any set of bands. 



5. Conclusions 

In this paper we have presented an optimized technique 
to produce highly accurate extinction maps from multi- 
band near-infrared photometric data. The method, which 



is a natural generalization of the near infrared color ex- 
cess method of Lada et al. (1994 ), is able to produce sig- 
nificantly less noisy (and thus more accurate) extinction 
maps taking advantage of all bands available. A first ex- 
ample of application of this new technique to 2MASS data 
has shown an improvement with respect to the standard 
Nice algorithm of a factor 2 on the noise variance. This 
way, we have been able to detect extended diffuse halos 
down to Ay ~ 0.5 magnitudes. 

Acknowledgements. We thank Charles Lada and Lindsay King 
for helpful discussions and observations on this work. We 
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Fig. 7. The same region shown in Fig. [5] but with the extinction map obtained using the standard method ( |Lada 
et al. 1994). A quantitative analysis of the noise, carried out in regions with the lowest column density, shows that the 
variance of this map is about a factor two larger than the one shown in Fig. |5|. Note also that, although the contours 
levels are the same as in Fig. ||, they appear much less smooth. This is another indication of the higher noise level of 
this map (rather than of intrinsic cloud structure; note that the resolution is the same for all images). 



of data products from the Two Micron All Sky Survey, 
which is a joint project of the University of Massachusetts 
and the Infrared Processing and Analysis Center, funded by 
the National Aeronautics and Space Administration and the 
National Science Foundation. 
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